Projected Quasiparticle Theory for Molecular Electronic Structure 
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We derive and implement symmetry-projected Hartree-Fock-Bogoliubov (HFB) equations and 
apply them to the molecular electronic structure problem. All symmetries (particle number, spin, 
spatial, and complex conjugation) are deliberately broken and restored in a self-consistent variation- 
after-projection approach. We show that the resulting method yields a comprehensive black-box 
treatment of static correlations with effective one-electron (mean-field) computational cost. The 
ensuing wave function is of multireference character and permeates the entire Hilbert space of the 
problem. The energy expression is different from regular HFB theory but remains a functional of an 
independent quasiparticle density matrix. All reduced density matrices are expressible as an integra- 
tion of transition density matrices over a gauge grid. We present several proof-of-principle examples 
demonstrating the compelling power of projected quasiparticle theory for quantum chemistry. 
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I. INTRODUCTION 

Although it has been more than 80 years since the 
formulation of Schrodinger's equation, there is still no 
black-box computationally efficient treatment of strong 
correlations in electronic structure theory. Strong cor- 
relations, also known as static or non-dynamic correla- 
tions, appear from exact and near-degeneracies in the 
wave function that render the one-determinant Hartree- 
Fock (HF) picture qualitatively incorrect. By "compu- 
tationally efficient" , we mean an effective one-electron 
theory with mean-field computational cost. Of course, 
exact diagonalization of the molecular Hamiltonian in 
the Hilbert space of the problem (a method known as 
full configuration interaction or FCl) includes all strong 
correlations but it is impractical except for the smallest 
systems because of its combinatorial computational cost. 
Other high accuracy models involve at least 0{M^) com- 
putational effort, where M is the number of orbitals. The 
method presented in this paper accomplishes the goal 
of being an effective one-electron theory with mean-field 
computational cost, accurately describing strong correla- 
tions (and more) of finite systems in a black-box manner. 

We derive and implement a wave function method 
named Projected Quasiparticle Theory (PQT) based 
on symmetry-projected Hartrce-Fock-Bogoliubov (HFB) 
equations. All symmetries (particle number, spin, spa- 
tial, and complex conjugation) are here deliberately bro- 
ken and restored in a self-consistent approach. We start, 
that is, with a symmetry broken single Slater determi- 
nant, and then restore the symmetry via projection to 
obtain a multi-determinantal wave function which is vari- 
ationally superior. Projection methods on a deformed 
(i.e. broken symmetry) HFB state have been used in 
nuclear physics for many years Unlike those used 
in quantum chemistry,— these projections are based on 
the generator coordinate methodji In particular, number 
projected Hartree-Fock-Bogoliubov (PHFB) is perhaps 
most widely used because the attractive character of the 
nucleon-nucleon interaction usually leads to HFB solu- 



tions that have lower energy than HF. The work that 
we here present builds upon the number projection for- 
malism of Sheikh and RingJi"— We expand and develop 
their technique to include many molecular symmetries 
not previously considered. 

PHFB is a variational problem where a deformed 
quasiparticle determinant |$) is optimized in the pres- 
ence of projection operators P which can represent a col- 
lection of symmetries. Sheikh and Ring^ realized that the 
energy arising from this minimization (see below) can be 
expressed as a functional of the HFB regular density ma- 
trix and anomalous density matrix. Thus, this problem 
can be solved by diagonalization of an effective Hamil- 
tonian matrix. Here, we extend their work and consider 
projections onto eigenfunctions of the particle number, 
spin rotation (both and Sz), complex conjugation, and 
spatial symmetry operators. In essence, we here present 
a method where all symmetries in the molecular wave 
functioii^'^ may be deliberately broken and variationally 
restored. This includes both continuous (infinite dimen- 
sion) such as spin rotation and discrete group represen- 
tations such as spatial symmetry and complex conjuga- 
tion; the latter two had not been considered previously. 
The resulting mathematical problem of variation-after- 
projection (VAP) seems formidable at first glance, yet as 
we show below, it can be converted into essentially one 
of diagonalization of an effective one-quasiparticle Hamil- 
tonian with mean-field computational cost. Alternatives 
to the diagonalization approach that we pursue are dis- 
cussed in the nuclear physics literaturci^'^ 

The research presented here was motivated by re- 
cent work in our group. In a series of papersii"— , we 
have proposed Constrained-Pairing Mean-Field Theory 
(CPMFT), a method for dealing with strong correlations. 
By introducing a fictitious attractive pairing interaction 
between electrons in an active space, CPMFT accounts 
for strong correlations and properly dissociates molecules 
into fragments with correct spatial and spin symmetries. 
CPMFT is black-box but does not have a wave function; 
it is in practice a one-particle density matrix functional 
whose two-particle density matrix is not N-representable 
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but whose energy can be optimized via an effective HFB 
diagonalization problemJ^ HFB breaks particle number 
symmetry, so CPMFT requires a chemical potential for 
controlling electron number. 

In its singlet-paired version^ CPMFT is intimately 
connected with Unrestricted HF (UHF)jii Its Lagrange 
multipliers for constraining pairing interactions within 
an active space inspired the resolution of an old prob- 
lem: how to obtain the ROHF wave function within an 
UHF framework.— This is essentially a constrained vari- 
ation problem of spin-projection of an unrestricted de- 
terminant. We should note that in open-shell cases, spin 
symmetry breaking occurs spontaneously when using un- 
restricted orbitals. On the other hand, when dealing with 
singlet states, it does not matter if the restricted HF 
(RHF) wave function is stable or unstable; in both cases 
the VAP wave function originating from a UHF guess 
becomes multireference in nature. This is akin to "delib- 
erate" symmetry breaking followed by a restoration step. 
In other words, spontaneous symmetry breaking (due to 
so-called "HF instabilities" ) is not a necessary condition 
in VAP schemes. Curiously, CPMFT can be interpreted 
as deliberate number symmetry breaking (triggered by 
HFB with a fictitious attractive interaction), followed by 
a symmetry restoration step generated by its energy and 
two-particle density matrix definitionj^ir— 

The lessons learned with CPMFT and our desire to 
develop a wave function theory for strong correlations, 
steered us into projection schemes. Number projected 
HFB (also known as projected Bardeen-Cooper-Schrieffer 
or PBCS) is widely discussed in textbookaii^ and yields 
the Antisymmetrized Geminal Power (AGP) wave func- 
tion of quantum chemistry, a model that has been ex- 
tensively studied] ^^1^^" — As shown here, projected HFB 
can be used not only to optimize AGP wave functions 
with singlet geminals (singlet pairing) but also to opti- 
mize broken symmetry geminals which, to the best of our 
knowledge, has not been done before. When all symme- 
tries are broken, the resulting projected wave function is 
of multireference character, one that permeates the entire 
Hilbert space of the problem. By this, we mean that ev- 
ery Slater determinant built from the natural orbitals of 
our broken symmetry determinant is allowed to overlap 
with our projected wave function. Remarkably, the pro- 
jected wave function accomplishes this feat with a linear 
number of parameters in an implicit factorization scheme 
driven by the variational principle. 

Viewed in their natural orbital basis, closed-shell 
singlet-paired AGP wave functions belong to the senior- 
ity zero sector of the Hilbert space of the problem and 
suffer from known drawbacks Seniority is here defined 
as twice the number of broken electron pairs in a deter- 
minant (i.e., the number of singly occupied spatial or- 
bitals) . For more details about the seniority coiicept and 
its use in electronic structure theory, see Ref. [23. Spin 
projection on top of number projection forces the wave 
function to permeate to all seniority sectors. By delib- 
eratively breaking spin symmetry and forcing a and /3 



spatial orbitals to be different, all electron pairs are bro- 
ken and spin contamination is introduced. In our case, 
however, spin projection operators in the HFB context 
restore the correct quantum numbers for UHF (S'^) and 
Generalized HF (GHF) (5^ and Sz) orbitals and density 
matrices. In GHF the spin orbitals are linear combina- 
tions of a and /? spin, yielding a density matrix with all 
spin blocks populated (usually referred to as noncoUinear 
configuration) . Access to the triplet pairing channel is re- 
quired if spin symmetry is to be broken and the broken 
symmetry wave function has HFB character We note 
in passing that the spin projection method developed in 
this work provides a solution to the extended Hartree- 
Fock method of Lowdin and Mayer— and Goddard's GF 
method^, albeit accomplished simultaneously with num- 
ber projection. 

Both number and spin projection operators (or their 
generalization in the case of non-singlet spin) can be writ- 
ten as integrals over gauge angles in the generator coordi- 
nate approach that we employ jii^ The restoration of these 
symmetries is mathematically accomplished here via dis- 
cretization of the gauge integral over modest size grids. 
The corresponding symmetries are U(l) for number and 
SU(2) [homeomorphic to SO (3)] for triaxial noncoUinear 
spin-projection. 

Restoration of discrete symmetries like point group 
(space) and complex conjugation yield wave functions 
that are linear combinations of the generalized AGP wave 
functions discussed above. Complex conjugation is an an- 
tiunitary symmetrj*2i and its restoration requires special 
consideration, as discussed below. 

The effective one-electron PHFB Hamiltonian includes 
both Fock and pairing pieces that depend on transition 
density matrices that are defined over the gauge grid. 
These transition density matrices yield via gauge inte- 
gration a correlated and factorizable two-particle den- 
sity matrix. Computationally, the most expensive step 
of PHFB requires that the two-electron integrals be con- 
tracted with the transition density matrices at every grid 
point, so the computational scaling of this generalized 
form of PHFB is comparable to the mean- field cost of sin- 
gle reference HF times the number of grid points, which 
we show below is fairly insensitive to the size of the sys- 
tem. Depending on basis set, system size, and energy 
gap, this computational cosIj^ will be 0{PP) to 0{M^). 
The entire myriad of linear scaling tools (including alter- 
natives to diagonalization) is of course available to help 
further reduce this computational cost to 0{M)M- The 
integration over grid points is trivially parallelizable. 

Remarkably, the PHFB wave function presented in this 
work can be obtained with mean-field computational cost 
even though its multireference character permeates the 
entire FCI space. The fundamental reason behind this 
feat resides in a projected energy expression that is of 
HFB form and can be solved for by a standard HFB 
diagonalization. 

From a mathematical perspective, the results here 
obtained are supported by an implicit coherent state 
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representation of non-orthogonal wave functions out of 
which the component of desired symmetry can be pro- 
jected while simultaneously satisfying the variational 
principlcji— The class of multireference wave functions 
studied in this work, which can be written as a quasipar- 
ticle density matrix functional, have remained completely 
unexplored in electronic structure theory. The only pre- 
viously known density matrix functional that has an ac- 
companying wave function is Hartree-Fock. 

The remainder of this paper is organized as follows. 
We discuss the main equations of our method, deferring 
details to the Appendix. A number of benchmark results 
are then presented and discussed followed by concluding 
remarks. 

II. THEORY 

A. General Considerations 

Fundamentally, our multireference wave functions are 
simply the result of projection from a broken symmetry 
quasiparticle determinant. Wc will provide a few neces- 
sary details of quasiparticle theory here; the interested 
reader should consult one of the many textbooks avail- 
able on the subjectiii^ The quasiparticle creation opera- 
tors that we use are obtained by a Bogoliubov-De Gennes 
transformation of the standard electron creation and an- 
nihilation operators: 

P\ - T.^U,,a] + V,,a,). (1) 
A quasiparticle determinant can then be written as 

M/2 
i=l 

where M is the dimension of the single-particle basis and 
I) is the empty state. Note that this quasiparticle deter- 
minant dwells in Fock space rather than in Hilbert space. 
In the Hartree-Fock-Bogoliubov method, one variation- 
ally minimizes the expectation value of the Hamiltonian 
with respect to the coefficients U and V defining the 
quasiparticle orbitals. Given the coefficients U and V, 
one can form the Hermitian density matrix 

p = V* (3) 

and the antisymmetric anomalous density matrix 

/t = V*U"^ = -UV^ (4) 

We remind the reader that the unprojected HFB wave 
function converges to HF in cases where the two-body 
interaction is repulsive, as in electronic structure. 

It may be instructive to write the HFB determinant in 
a slightly different manner. In the so-called HFB "canon- 
ical" (natural orbital) basis, the quasiparticle determi- 
nant is 

s 

ii>)=AAn(i+^fc44.)i)' (5) 

fc=l 



where is a normalization factor and Cfe = Vk /uk ■ The 
product runs over the s orbitals defining the subspace 
over which the HFB wave function is allowed. Usually, 
s = M/2, but this is not always the case. Indeed, if 
s ~ N/2 where N is the number of electrons, then the 
number projected HFB corresponds exactly to HF. Note 
that s is known in the AGP literature as the rank of the 
geminal. In writing Eqn. [SJ we have used the indices k 
and k to represent the orbitals that are paired. Typically 
the paired orbitals are chosen to be ka and fc/3, the two 
spin orbitals formed from the same spatial orbital. When 
this is so, the projected HFB wave function is manifestly 
of seniority zero. However, one can choose to pair orbitals 
differently, a point to which we shall return later. 

From the form of Eqn. [S]it is clear that the HFB deter- 
minant contains contributions from states of many differ- 
ent particle numbers, and that for each particle number, 
number projection of HFB yields a linear combination of 
many determinants. In fact, when s = M/2 the number 
projected unrestricted HFB yields a linear combination of 
every determinant, which is what we mean when we say 
that the projected HFB wave function permeates Hilbert 
space. Pairing the spatial orbitals ka and kf3 instead 
yields a projected HFB wave function which is a linear 
combination of every determinant of seniority zero. In 
either case, each determinant has its own coefficient, but 
those coefficients are determined by the HFB parame- 
ters Cfe. In other words the coefficients of the different 
determinants are factorizcd. 

Approximate wave functions such as quasiparticle de- 
terminants need not possess the same symmetries as the 
exact solutions. In fact, constraining the variationally 
optimized state to preserve certain symmetries can only 
lead to higher energy solutions as this reduces the varia- 
tional manifold, a fact now known as Lowdin's symmetry 
dilemma^ 

Projection operators can be used to restore symme- 
tries of the Hamiltonian from a broken-symmetry ap- 
proximate wave function, thus avoiding the symmetry 
dilemma. One can use the projection operators in two 
ways: 

• In the projection-after-variation (PAV) scheme, one 
optimizes a broken symmetry wave function, and 
then performs a single-shot symmetry restoration 
with the projection operator. 

• In the variation-aftcr-projection (VAP) scheme, the 
optimization is carried out in the presence of the 
projection operator. That is, one minimizes the ex- 
pectation value of the projected wave function with 
respect to variations of the underlying deformed 
state. 

In most cases VAP is preferred since it uses all the varia- 
tional flexibility available and avoids the artifactual phase 
transitions associated with symmetry breaking of the un- 
derlying wave function, which may become more pro- 
nounced when using the PAV schemei^i^ Nonetheless, 
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carrying out the full variational optimization with the 
VAP scheme may lead to equations that become signifi- 
cantly more involved than solving the original variational 
problemj^i^S In finite systems, quantum fluctuations in- 
corporated in the VAP scheme lead to multireference 
wave functions that remove all artifactual phase tran- 
sitions due to spontaneous symmetry breaking. This will 
be clear in the dissociation curves presented below. We 
note in passing that in the thermodynamic limit (infinite 
systems), VAP and PAV based on a mean- field state do 
not improve the energy over that of the broken symme- 
try mean-field solutioui^i^ The focus of this paper is on 
finite systems only. 

In the special case when the deformed state is a quasi- 
particle determinant, it turns out that the VAP equa- 
tions, though involved, can be solved in a straightfor- 
ward manner. The nuclear physics literature contains 
a large number of approximate solutions^ Exact solu- 
tions to certain specific VAP problems in nuclear physics 
have also been discussed^i^ii^ but not always fully im- 
plemented. The projection operators employed here are 
based on the generator coordinate methodii^"— and are 
very different from the traditional approach in quantum 
chemistry due to L6wdin<^ We note that the "Variation 
After Mean field Projection In Realistic model spaces" 
of Schmid and collaborators share common features with 
out formalism, though the authors do not optimize the 
projected HFB state via diagonalization of a Hermitian 
matrix as we do.— To the best of our knowledge, projec- 
tion methods have been applied to the electronic problem 
only twicCf^i^ and then only with model Hamiltonians. 

Given a Hermitian operator A which commutes with 
the Hamiltonian, we can choose eigenfunctions of an ap- 
proximate Hamiltonian such as HF or HFB to also be 
eigenfunctions of A. If we break A-symmetry (that is, if 
our wave function is chosen to not be an eigcnfunction 
of A) , we can construct a corresponding unitary operator 
U = exp i0(A — A) with parameters and A where (jj is 
a gauge angle over which we will integrate and A is the 
desired eigenvalue of A. Integration over will then yield 
a projection operator P which projects eigenstates of A 
with eigenvalue A out of the the symmetry-broken wave 
function, provided that all eigenvalues of A are rational. 
We sketch the main ideas below for particle number pro- 
jection as an example. 



B. Particle Number Projection as an Example 

Suppose that |$) is a quasiparticle determinant. As a 
quasiparticle determinant, it is completely specified by its 
density matrix p and anomalous density matrix k, and 
is a linear combination of particle number eigenstates: 



where N = ^ajai is the number operator. We empha- 
size that the number eigenstates are multidetermi- 
nantal wave functions rather than being single determi- 
nants. 

Then consider the unitary operator l}{9) — cxp(i6'iV). 
Its action on |<I>) is to produce other quasiparticle deter- 
minants: 

C/(0)|$)=^c..e'«^'^-|*,)^|$(0)). (7) 

We can project out the component \^ j) with parti- 
cle number Nj by multiplying by the weight function 
Wj{0) = ex-p{—i6Nj)/{2TT) and integrating. That is, 

c,|*,)= J d0w,ie)Ui9m. (8) 



We will write 

27r 

Pj^ I dew,{9)Ui9), (9) 



which is the projection operator onto the number eigen- 
state of interest. 

Now, the energy we wish to minimize is simply 

mp]p,\^) imi^) 

where we have used the fact that Pj is Hermitian, idem- 
potent, and commutes with the Hamiltonian. The matrix 
elements in the numerator and denominator can be, in 
principle, expressed entirely in terms of the projection 
operator and the density matrices p and k associated 
with 1$). We then simply minimize E with respect to 
p and K subject to the constraint that they arc com- 
patible with a quasiparticle determinant, and the result 
defines the VAP energy and wave function. Projections 
onto eigenstates of other operators follow a qualitatively 
similar pattern, though with significant differences in the 
details. 

In general, given a Hermitian constant of motion (sym- 
metry) , it is always possible to construct a unitary opera- 
tor like that in Eqn. [T] This operator creates a manifold 
of degenerate states from which the desired state of in- 
terest can be extracted via a projection operator of the 
form in Eqn. ^ It is not the purpose of this paper to dis- 
cuss these mathematical tools in detail. The interested 
reader is referred to one of the textbooks in the field lii^ 



C. Projection Operators 



1$) = Cfel'I'fe), (6a) Here we wish to briefly review the projection operators 

we will use. We wish to take a fairly general form, to 
A^l^/s) — A^fel^'fe), (6b) make what follows as transparent as possible. 
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When the generators in which we are interested are 
continuous, the projection operators we consider can be 
expressed as 



P,(J) 



de* 



(11) 



where we wish to project onto the j"^ cigenstate of the 
operator J. The precise form of the weight j{9) de- 
pends on the operator J and also on the eigenstate j. 
The rotation operator R{9, J) similarly depends on the 
symmetry operator J. If the generator is instead discrete, 
the integration is replaced by a summation. 

The rotation operator R{9, J) acts on a deformed or 
broken symmetry state in such a way that 



|($|$)|= {<^{e,jme,j)) 

(ci>|i|$) = ($(0,J)|i|$(0,J)) 



(12a) 
(12b) 



where |<I>(0, J)) = R{9,J)\(^) and where A commutes 
with J and thus with R{9,J). In other words, R{9,J) 
acts on the deformed state, preserving the norm and ma- 
trix elements of commuting observables up to an over- 
all phase factor. Typically, R{9, J) will be unitary, 
but it may instead be antiunitary. We point out that 
the rotated states form a nonorthogonal set that can 
be overcomplete, in the same manner as are coherent 
states. In fact, as we have already noted, the success of 
this approach is based on an underlying coherent state 
representation^^ 

In this paper we will consider four types of projection 
operators associated with particle number, spin rotation, 
point group symmetry, and complex conjugation. We 
will briefly discuss the form of each of these projection 
operators. We will also give the matrix representations 
R(0, J) of the rotation operators R{9, J) in an orthonor- 
mal basis of spin orbitals \ia), where we have written 
them in spin blocks. That is. 



R,,i9,J)^mRi9j)\j). 



(13) 



Number - The projection operator associated with 
particle-number restoration is particularly simple. It is 
given as an integration over the gauge angle (j), in the 
form of Eqn. [Til with 



WN{(f>) = 


1 -i<pN 

2tt 


(14) 


Ri<P,N) = 




(15) 


R(0,iV) = 


[ e'-^li 


(16) 



where N is the number operator and N is the number of 
particles that the projected wave function shall possess. 

Spin - The electronic Hamiltonian is spin free, and 
thus the wave function should be invariant to rotations 
of spin {i.e. spin can be quantized along an arbitrary 



axis). Formally, this procedure corresponds to using a 
projection operator of the form 



M.K 



(17) 



where S refers to the eigenvalue of S*-^, M to the eigen- 
value of Sz, and the cm are variationally optimized co- 
efficients. Because we do not generally know the spin 
eigenfunctions \S;M), we instead write the projector as 
an integration over spin rotations characterized by the 
Euler angles il = (a,/3,7). That is, we have 



Ps 



CMC*K / dnwsM.KmR{n,S) 



(18) 



where the weights and the rotation operator are given by 



2S + l ^s. 
87r2 



(19) 
(20) 



R{n,S) = e'"'^' e'^'^y 
R(n, S) = R(a, ^.) R(/3, Sy) R(7, 5,), (21) 



R(a,i>^) = ( ^ ^_i„/2 ^ 



c-'" 
cos(^/2) 1 sin(/3/2) 1 



(22) 



R(/3, Sy) - sin(^/2) 1 cos(^/2) 1 ) ' ^^^^ 

Here, Dfjj^{n) = {S; M\R{n, S)\S; K) are Wigner rota- 
tion matrices. "44. 

We note that the spin projection operator has the 
same form as the angular momentum projection oper- 
ators commonly used in nuclear physics to restore spa- 
tial rotations, as the SU(2) algebra characterizing spin is 
homcomorphic to the SO (3) algebra characterizing angu- 
lar momentum. We refer the reader to the exposition of 
the angular momentum projection operators in Ref. [ll. 
Note that if the underlying reference state is an eigen- 
function of Sz (a collinear state), as is usually the case, 
then the integrations over a and 7 are trivial and the 
projection operator becomes 



^S,Ns,Ns — — ; 



where Ns 



1 



'd/3sin(/3)4^^^(/3)e''^^^ (24) 



and d%^j^^{/3) = 

{S;Ns\R{l3,Sy)\S;Ns) is Wigner's smah d-matrix4i 
Note also that in restoring invariance to the axis of spin 
quantization we also force the wave function to be an 
cigenfunction of S^. This form of the spin projection 
operator was proposed by Percus and Rotenberg^ 
in 1962 and was used in chemistry by Lefebvre and 
Prat^ in the context of HE, though their methods and 
techniques never became widely used. 

Point Group - Restoring point group symmetry can 
follow the same mathematics as used above but is more 
simply accomplished by diagonalization of the Hamil- 
tonian matrix in the basis of the wave functions /t |$), 
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where /i are the elements of the point group under con- 
sideration. In other words, the projected energy becomes 



(25) 



where are variational coefficients. This is similar to the 
way other authors restore parity in nucleiiiS The rotation 
matrices for spatial symmetry restoration correspond to 
the matrix representation of the symmetry elements of 
the point group. 

Spatial symmetry breaking and restoration is useful for 
exploring regions of potential energy surfaces where some 
spatial symmetry is preserved, and only those symmetry 
operations which are preserved should be restored in the 
manner described above. If different symmetry operators 
are restored at different nuclear geometries, the result- 
ing potential energy surface is likely to be discontinuous. 
One could try to extend our methods by using all the 
symmetry operators defined at high symmetry points in 
the diagonalization above, even at points of lower sym- 
metry, but this is not symmetry breaking and restoration 
and its discussion is beyond the scope of this paper. 

Complex Conjugation ~ The exact wave function in 
quantum chemistry can be chosen to be an eigenfunction 
of the complex conjugation operator K as long as the 
Hamiltonian commutes with it^ which is always true in 
the nonrelativistic case where H is real. Complex conju- 
gation is an antiunitary operator so its spectrum does not 
carry good quantum numbers with iti^ A wave function 
l^") = J2 is an eigenfunction of K when 



(26) 



D. Variation- After-Projection Scheme 

In the variation-after-projection scheme, our aim shall 
be to minimize the expectation value 

^ {mP.m ^ jd9w,i0){mRiom ^ ^28) 



Jd0w,{e){<i>\R{0m 



The minimization should be carried over all possible 
states 1$) of the HFB form. 

It is most straightforward to evaluate Hamiltonian el- 
ements between quasiparticle determinants in intermedi- 
ate normalization. We can define such determinants as 



1^) 



mm 



Then the energy expression is 

jdew^iemiimm 



jd9wj{e){<!>\R{e)m 

Jdew.je) imiem {m\o) 
jdew,{e)mRi9m 



(29) 

(30a) 
(30b) 



Following Sheikh and Ring^, we make this expression 
more compact by absorbing the overlap matrix elements 
into new weight factors, defining 



in terms of which the energy is simply 

E,= f dey.iemm. 



(31) 



(32) 



implying 



(27) 



for some real number cf). An HFB wave function need not 
be an eigenfunction of K, and if the HFB wave function 
does not satisfy the property of Eqn. [^H then we can di- 
agonalizc the Hamiltonian in the basis {l^*}, The 
resulting wave function has lower energy. 

A cautionary note is in order. The matrix elements in- 
volved in the CI problem required to restore point group 
symmetry can be found by building the rotation matrix 
corresponding to the operator jj, and using the formulas 
provided below. On the other hand, the matrix elements 
involving the complex conjugation operator are slightly 
more involved as there is no associated rotation operator. 

Henceforth, we will suppress dependence on J for 
brevity of notation. 



E. Evaluation of Matrix Elements 

The remaining task is to evaluate the Hamiltonian and 
overlap elements. The overlap matrix elements can be 
expressed asii^ 



imim) = exp { ^Tr [log (1 - M{9))] 
exp<j iTr[log(l-M(0))] 



(33) 



Here and in what follows, we use the matrix 

M{0) = Z*R(6')ZR"^(6I), (34) 

where Z = V*(U*)~^ is the Thouless matrijiii^ charac- 
terizing the HFB state |4>) in terms of its orbital coeffi- 
cient matrices U and V. We point out here that we have 
considered only single-particle operators, which places re- 
strictions on the symmetries to be projected. Note also 
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that we will frequently write in place of R(0) in or- 
der to reduce clutter, and so on for other ^-dependent 
matrices. 

The Hamiltonian matrix elements (^\H\9) = Ha are 



+ \Y.^ij\\ki)R*^{e)Kki{e). 



(35) 



Here, hij are the usual one-electron integrals and {ii\\kl) 
are antisymmetrized two-electron integrals in Dirac nota- 
tion, while pg, Kg, and Kg are transition density matrices 
given by 



{^\a,a,\e) 
(Z*)-iMe(l- Ale)"' 



-{l-Me 



(36a) 
(36b) 
(36c) 
(36d) 

(36e) 
(36f) 



One can simplify the Hamiltonian matrix elements as 



He ^ ^Ti [{h + J^g)pg - A*gKe] 
= ^TT[{h + Te)pg-AeK*g], 



(37) 



where J^g is a generalized Fock operator and A.g and 
A.g are generalizations of the pairing matrix of HFB and 
CPMFT: 

T^kiO) = h,,, +J2Mkl) PI J (0) (38a) 
= h,k+G,kiO) (38b) 

^*.^(^) = ^E<*^1l^^)'*fe'(^) (38c) 

^^^■(^) = ^E<*-?1l^^)^fc'(^)- (38d) 

The overlap matrix elements and the transition den- 
sity matrices can be explicitly written in terms of the 
density matrix p and the anomalous density matrix k of 
the underlying HFB state |(f>)»i In particular, the overlap 
matrix elements are written as 

where the matrix Ce is constructed from p and k accord- 
ing to 

Cg 1 = pRgpRl - kRIk*rI. (40) 
The transition density matrices are in turn expressed as 

pg = RepRjCep, (41a) 

(41b) 
(41c) 



Kg = RgpRlCgK, 
K*g = R*gK*RlCgp. 



Since the matrix elements are all functionals of p and k, 
so too is the projected HFB energy itself: 



Ej =Ej[p,k] 



(42) 



We remind the reader that the index j labels the quantum 
numbers which have been projcctively restored. 



F. Variational Equations 

The importance of the foregoing result cannot be 
overemphasized. The fact that all matrix elements can 
be expressed in terms of the density and pairing matri- 
ces of the underlying HFB state implies that one can 
minimize the energy directly with respect to p and k, 
obtaining an effective mean-field Hamiltonian analogous 
to the Fock operator or the quasiparticle Hamiltonian of 
HFB. Specifically, one can optimize the functional 



C[p, k] = Ej [p, k] - Tr [A (7?.^ - 71}] 



(43) 



where A is a matrix of Lagrange multipliers used to con- 
strain the generalized density matrix TZ-, expressed in the 
Valatin form as 



n 



p K 

-K* 1 p* 



(44) 



to remain idcmpotcnt. The idempotency of the gener- 
alized density matrix is equivalent to the requirement 
that the state remains a quasiparticle Slater determinant. 
Making the functional stationary with respect to varia- 
tions in 7?. leads to the condition 



[^,,7?.] =0, 



(45) 



which is simply the Brillouin condition for HFB. Here, 
the effective Hamiltonian matrix Hj is given by 



SE^ 
STZ' 



(46) 



Expressions for the matrix elements of Hj are provided 
in appendix A. 

Equation l45l thus constitutes the form of the projected 
HFB equations expressed in matrix form. Solving the 
PHFB equations is thus conceptually simple. Given an 
initial guess of the quasiparticle density matrix 7?., one 
constructs the effective Hamiltonian matrix fij, diago- 
nalizes it, and constructs an updated IZ using the eigen- 
vectors of Tij. Convergence is achieved once Eqn. 1151 is 
satisfied up to a previously determined tolerance. 

There are some subtle differences between the pro- 
jected and the regular HFB variational problems that 
merit a few words. In the regular HFB equations, one 
must introduce the chemical potential /x forcing the HFB 
determinant to contain the correct particle number on 
average. This chemical potential is not required in pro- 
jected HFB, but we include it nonetheless as it often im- 
proves convergence of the projected HFB equations. The 
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chemical potential, to be clear, is applied to force the 
quasiparticle determinant to contain the correct number 
of electrons, which improves convergence of the projected 
HFB equations but does not change the final result. 

There is also the key question of how one should occupy 
the quasiparticle orbitals after diagonalization of the pro- 
jected HFB Hamiltonian. This has been discussed in the 
literature in the case of HFB 4^ We emphasize that oc- 
cupying the lowest-energy orbitals in the projected HFB 
Hamiltonian spectrum need not lead to convergence, or 
to the lowest energy solution even when convergence is 
achieved. Using the fact that, if the initial guess is good 
enough, then Eqn. |45]is approximately satisfied at every 
cycle, we diagonalize the modified Hamiltonian matrix 



(47) 



with A < 0. When the equations are converged, the 
spectrum of fij is identical to that of Hj except that the 
occupied quasiparticle orbital energies have been shifted 
down by A. This "level shifting"— allows us to occupy 
the desired orbitals following an aufbau principle. At 
self-consistency, it is of course irrelevant whether we have 
forced to commute with H or with H. 



G. Simplified Forms of tfie Variational Equations 

We have written the projected HFB equations, Eqn. |45] 
in the spin-orbital framework, allowing for quite general 
symmetry breaking. However, the HFB and PHFB equa- 
tions, like the HF equations, have self-consistent sym- 
metries. In other words, if we choose an initial den- 
sity matrix to contain certain symmetries, then the self- 
consistent density matrix will also contain those symme- 
tries. Thus, for example, if we were to supply projected 
HFB with an RHF initial guess, we would converge to 
an RHF solution as RHF is a special case of the more 
general projected HFB. It is up to the user, therefore, 
to decide which symmetries, if any, are to be broken in 
the unprojected HFB state. Once this choice has been 
made, the procedure for restoring these symmetries is 
fairly general. 

We can impose self-consistent symmetries on the un- 
derlying HFB state to simplify the PHFB equations, and 
discuss two such forms which we have used in this paper. 
Yamaki et al^ have discussed other possible structures 
of simplified problems in the context of regular HFB. 

In the restricted HFB (RHFB) equations with singlet 
pairing, we prepare our initial guess for the quasipar- 
ticle density matrix 72. imposing the following condi- 
tions: paa = Pl3l3, Pa/S = Pfja = 0, Kaa = K^fSP = 0, 

'*a/3 = {i^ap)~^ ■ Note that this form is appropriate for 
closed shell systems. One can then solve the projected 
HFB equations using only half the dimension with the 
simplified quasiparticle density matrix 7?., given by 



and the simplified effective Hamiltonian 

TV j — — . 



(49) 



Generally the RHFB wave function does not break spatial 
symmetry, but it retains the flexibility to do so. 

The unrestricted HFB (UHFB) equations allow for 
spin symmetry breaking, though the overall HFB state 
remains an eigenfunction of Sz- Here, we prepare our ini- 
tial guess for the quasiparticle density matrix assuming 
the following conditions: Pa/s = ppa = 0, Kaa = i^pi3 = 
0. This leads to the following structure for the simplified 
quasiparticle density matrix: 



(50) 



It = 



paa. 



(48) 



This allows for opposite spin (7715 = 0) triplet pairingi^ 
in the underlying quasiparticle determinant while ex- 
cluding same spin (mg = ±1) triplet pairing. Break- 
ing the Sz symmetry of the underlying HFB state intro- 
duces same-spin triplet pair correlations to the underly- 
ing quasiparticle determinant, and requires solving the 
PHFB equations in the general framework, which we de- 
note as GHFB. 



III. COMPUTATIONAL DETAILS 

We have implemented the projected Hartree-Fock- 
Bogoliubov equations in the way described above both 
in an in-house code and as part of the GAUSSIAN^ suite 
of programs. We have validated our implementation of 
the particle-number projected equations by comparing 
results with an AGP code that we had available in our 
research group.— We compared the energy, the natu- 
ral orbital occupations and the geminal coefficients and 
reached quantitative agreement in all cases tested. 

In this work, we only report tests on a few small sin- 
glet systems, using basis sets of minimal to double-^ + 
polarization quality. Calculations on larger systems and 
with larger basis sets are possible, but we present here 
pilot calculations to illustrate the main features of the 
method. 

We have a wide variety of possible states, depending 
on which symmetries we wish to break in the underly- 
ing quasiparticle determinant and which broken symme- 
tries we wish to restore by projection. We summarize 
our nomenclature in Table HI Note that we have not 
exhausted all possible structures of the quasiparticle de- 
terminant |<I>) and. therefore, the density matrices p and 
1^50^ The versions we have used closely resemble the forms 
of Hartree-Fock determinants most commonly employed. 

Let us pause to discuss the hierarchy of wave functions 
that fall into our classification. Starting from a broken 
symmetry HFB determinant, number projection results 
in a single AGP wave function, with broken symmetry 
if the quasiparticle determinant is of UHFB or GHFB 
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TABLE I. Nomenclature defining our various projected and 
unprojected HFB states. 



Designation Form of p 



Form of k 



RHFB 

UHFB 
GHFB 



paa 

Paa 
' paa 

Pan 

Polol pap 

yPl3cy PflfS 






l^afj 




a 









) 


h!>aa 


l^aP \ 


K-I3a 





Designation 


Symmetry Restored 


N 


Particle Number 


S 


Spin 


K 


Complex Conjugation 


Ci 


Inversion Point Group Symmetry 



character. Adding spin projection yields a linear combi- 
nation of AGP wave functions such that the overall result 
is a spin eigenfunction, and restoring discrete symmetries 
then yields a linear combination of these wave functions. 
It is not yet entirely clear how one could a priori deter- 
mine which symmetries must be broken and restored in 
any particular calculation or which of these symmetries 
is more important in one molecule than another. Future 
calculations will surely clarify this and other aspects of 
the methodology that we present here. 

We prepare the initial guess to our restricted number- 
projected (NRHFB) calculations by diagonalizing the 
core Hamiltonian (neglecting the electron-electron repul- 
sion) or by solving restricted HF equations. We use 
the orbital energies obtained to thermalize our initial 
guess of occupations according to a Fermi-Dirac distri- 
bution. Wc then build the singlet-pairing Ka/3 matrix 

using Kai3 = Paa " Pqq- 

For our spin symmetry broken calculations, our cur- 
rent approach is simple: we converge the restricted equa- 
tions, and then mix the highest-occupied and the lowest- 
unoccupied quasiparticle eigenvectors using some pre- 
defined angle. Constructing the quasiparticle density 
matrix from such a set of eigenvectors leads to a spin- 
symmetry broken guess of 7?., that we then iterate to 
convergence. The quality of this initial guess is rather 
poor, and improved initial guesses should help alleviate 
the convergence difficulties described below. 

The number of grid points required in the discretiza- 
tion of the particle-number and spin projection inte- 
grations is relatively low. When using the simplified 
form of the equations appropriate for particle-number 
projection)^ we have observed that convergence of the 
energy to 1 nano-Hartree is achieved with 7 or 9 grid 
points. Use of the more general equations presented in 
this paper seems to require a larger grid even for particle 
number projection, and we have used 15 points for most 
of our calculations. Discretization over the Euler angles 
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FIG. 1. Dissociation of LiH in the cc-pVDZ basis. NRHFB is 
equivalent to AGP, and offers sizable improvement over RHF 
and UHF. 



in spin projection leads to good convergence of our cal- 
culations when using about 10 grid points per angle. We 
generally integrate with equally spaced grid points using 
the trapezoid rule but use Gauss-Legendre quadrature 
for the integration over /? in the triaxial spin projection. 
We discuss the accuracy of the integration grid in Section 
IIV Gl We emphasize that the calculations at each grid 
point are independent, and trivially parallel. 

Finally, we should point out that the self-consistent 
equations of PHFB can be rather difficult to converge 
and we have employed the direct inversion of the iterative 
subspace (DIIS) algorithm^i^ in addition to the level 
shifting described previously. 



IV. RESULTS AND DISCUSSION 

A. Number Projection and AGP 

The simplest case of projected quasiparticle theory 
is number projection of a restricted HFB determinant, 
which we refer to as NRHFB and which is identical to 
AGP with singlet geminals. As mentioned earlier, we 
have used this fact to test the correctness of our NRHFB 
implcmention. More interestingly, however, wc can use 
our NRHFB code to generate AGP wave functions at 
mean field cost. To the best of our knowledge, the 
previously most efficient formulations of AGP scale as 

We note that the NRHFB/ AGP wave function dwells 
strictly in the seniority sector of Hilbert space and is a 
linear combination of each and every seniority determi- 
nant with, however, restrictions relating the coefficients. 
In other words, every determinant in the NRHFB is a 
closed shell, and each closed shell determinant is assigned 
a coefficient. 
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To illustrate the efficacy of NRHFB, wc show the dis- 
sociation of LiH in Fig. [TJ As it is well known, AGP is 
exact for one electron pair, so LiH is perhaps the sim- 
plest dissociation case for which the method is not FCI. 
As always, the RHF wave function is unable to dissociate 
to open-shell fragments, and the broken symmetry, spin- 
contaminated UHF solution splits off from it and disso- 
ciates to two UHF atoms. Meanwhile, NRHFB offers 
significant dynamic correlation at equilibrium and disso- 
ciates cleanly to two ROHF atoms. The CASSCF(2,2) 
is the minimal CAS needed to dissociate the molecule 
properly to ROHF fragments, and is included for com- 
parison. In this particular case, NRHFB is below the 
minimal CAS, although as we shall see this is not always 
so. 



B. Unrestricted HFB and AGP 

Perhaps the simplest extension to NRHFB is to allow 
the breaking of spin symmetry to give us NUHFB. This is 
equivalent to an AGP with unrestricted (and potentially 
broken symmetry) orbitals. As pointed out by Weiner 
et aZj^, this is likely to be beneficial for the description 
of molecular dissociation, but to the best of our knowl- 
edge this is the first time that variationally optimized 
AGP with broken symmetry orbitals has been reported. 
We should also add that simply by breaking spin sym- 
metry, we allow the NUHFB wave function to permeate 
all of Hilbert space. In other words, it spreads beyond 
the seniority sector and gives us contributions from, in 
principle, every Slater determinant with the appropriate 
numbers of a-spin and /3-spin electrons. 

In Fig. [5] we show the dissociation of C2H4 to two 
triplet CH2 fragments. Because we are breaking a dou- 
ble bond, one can expect correlation effects to be very 
important. One consequence is that the UHF curve sep- 
arates from the RHF curve near equilibrium. A second 
consequence is that NRHFB ends up between the RHF 
and UHF limits (though much closer to the latter). In 
order to get a reasonable description of the dissociation, 
we arc forced to break spin symmetry to give us NUHFB, 
which goes to a dissociation limit below two ROHF frag- 
ments. Interestingly, the NUHFB curve is nearly parallel 
to the UHF curve, except near equilibrium where it picks 
up a little extra correlation. We suspect that the dissoci- 
ation limit amounts to one UHF fragment and one AGP 
fragment, though we cannot confirm this as we do not at 
present have the capability to do open-shell AGP calcu- 
lations. The minimal CAS is also shown; in this case it 
provides a slight variational improvement at equilibrium 
and dissociates to two ROHF fragments. We can mea- 
sure the spin contamination introduced by the broken 
symmetry NUHFB wave function, and Fig. [2] also shows 
the expectation value of 5'^ using the formulas provided in 
Appendix B. It is interesting to note that for large bond 
lengths the NUHFB and UHF wave functions contain 
identical amounts of spin contamination. The bumps in 
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FIG. 2. Top Panel: Dissociation of C2H4 to two CH2 frag- 
ments, in the cc-pVDZ basis. Symmetry-broken AGP {i.e. 
NUHFB) offers significant variational improvements. Bottom 
Panel: Spin contamination in UHF and NUHFB as a function 
of C-G bond length. 
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FIG. 3. Torsion of C2H4 about the C-G double bond, in the 
cc-pVDZ basis. Breaking spatial and spin symmetry offers 
significant amounts of static correlation at both the HP and 
projected HFB levels. 



the UHF and NUHFB curves are manifestations of spa- 
tial symmetries being broken as the bond is stretched. 
We note that the physical dissociation of C2H4 to triplet 
CH2 fragments is accompanied by substantial changes in 
the geometry of the CH2 fragments which wc have here 
ignored but which are essential if one is to accurately 
reproduce experimental results.— 

A second example where spin symmetry breaking plays 
an important role is in the torsional barrier of C2H4, 
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shown in Fig. [3] This is a familiar muhireference prob- 
lem, where the correct curve requires two determinants at 
= 90°. The RHF curve we have shown predicts a very 
large barrier, which is greatly reduced by UHF. Similarly, 
the barrier in NRHFB is predicted to be rather large, but 
much smaller with NUHFB which again closely parallels 
the UHF curve. The minimal CAS predicts the barrier 
to be somewhere between that predicted by NRHFB and 
NUHFB. Note that we have ignored non-adiabatic effects 
which should be relevant here. 

In order to provide a more formal description of what 
we mean by breaking spatial and spin symmetry in the 
context of HFB, we return to the form of the quasipar- 
ticle determinant, as written in Eqn. [5l We recall that 
the orbitals to be paired, represented by indices k and 
k, are usually the a and /? spin orbitals corresponding to 
the same spatial orbital. Breaking spin symmetry means 
relaxing this constraint. In particular, when we consider 
general spin orbitals^ to construct the quasiparticle vac- 
uum, the spin label is lost as the orbitals have contribu- 
tions from both a and j3 character. In such a case, it 
is still true that rifc = n^., and those correspond to the 
orbitals paired, but it is no longer true that one of the 
paired orbitals is a and the other /?. In breaking the 
spin symmetry of the wave function we have broken the 
time-reversal character of the Cooper pairs. 



C. Restoring Spin Symmetry 

Having broken spin symmetry, the next step is to re- 
store it, yielding what we caU SNUHFB and SNGHFB. In 
both cases the projected HFB wave function is an cigcn- 
function of S"^ and of Sz- The distinction between them 
is that the unprojcctcd UHFB wave function is an eigen- 
function of Sz, while the unprojected GHFB wave func- 
tion is not. We emphasize that the RHFB wave function 
is an eigenfunction of both 5^ and of Sz] consequently, 
spin projection on the RHFB wave function merely re- 
turns the RHFB wave function. 

Figure |3] shows the dissociation of N2. As is well 
known, the RHF curve is disastrous, while UHF strongly 
underbinds and predicts a barrier to the formation of the 
N-N bond. As with the dissociation of C2H4, NRHFB 
goes to a dissociation limit between the RHF and UHF 
limits while offering some improvement at equilibrium. 
In this case, NUHFB follows NRHFB near equilibrium 
but goes to the UHF limit at dissociation; that NUHFB 
does not go below the UHF limit at dissociation is sim- 
ply because we are working in a minimal basis. Restoring 
the broken spin symmetry in SNUHFB offers substantial 
improvements all across the dissociation curve while still 
dissociating to two ROHF atoms. 
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FIG. 4. Dissociation of N2 in the STO-3G basis. Restoring 
the broken spin symmetry provides significant correlation and 
eliminates the instability in the NRHFB wave function. 



D. Restoring Complex Conjugation Symmetry 

Figure 2] also shows the dissociation curves predicted 
by KNRHFB and KNUHFB for the nitrogen molecule. 
Complex conjugation restoration recovers a very sig- 
nificant fraction of correlation near equilibrium. KN- 
RHFB does not dessociate to the correct limit, a fea- 
ture that KNUHFB solves by breaking spin symmetry. 
Adding complex conjugation on top of spin projection 
(KSNUHFB) yields a curve that stays very close to full 
CI for all bond lengths. 

A particularly interesting test case is provided by an 
equally spaced linear chain of hydrogen atoms. We here 
limit ourselves to H4 in a minimal basis. As can be seen in 
Fig. [5l RHF dissociates incorrectly while UHF yields four 
hydrogen atoms. As we have come to expect, NRHFB 
(even with broken spatial symmetry) gives a dissociation 
limit between RHF and UHF, while NUHFB goes to the 
UHF dissociation limit. 

Restoring additional symmetries leads to results nearly 
identical to full CI, and Fig. |6] therefore shows the 
deviation from full CI. We achieve significant improve- 
ment over NUHFB by restoring coUinear spin symme- 
try. Restoring noncollinear spin symmetry in NGHFB 
appears to yield results agreeing with full CI to better 
than the nano-Hartree level, though we have experienced 
great difficulties in converging the equations and have 
been unable to generate the entire potential energy curve. 
We speculate that the accuracy of the SNGHFB is tied in 
some way to the fact that we have only two symmetry- 
unique hydrogen atoms (and thus only two symmetry- 
unique electrons) and are working in a minimal basis. 

Rather than restoring spin symmetry, we can choose 
to restore complex conjugation symmetry, and in so do- 
ing achieve the correct dissociation limit with restricted 
orbitals. 
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FIG. 5. Dissociation of equally spaced linear H4 in the STO- 
3G basis. Only by breaking symmetry can projected HFB 
reach the correct dissociation limit. 
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FIG. 6. Comparison of KNHFB and SNHFB to full CI for the 
dissociation of equally spaced linear H4. Restoring these sym- 
metries leads to results almost equal to full CI. In particular, 
SNGHFB appears to be energetically identical to full CI for 
the points at which we successfully converged the equations. 



E. Restoring Spatial Symmetry 

Throughout our discussion, the reader may have no- 
ticed that our projected HFB wave functions generally 
dissociate to ROHF fragments. Of course projected HFB 
calculations on the isolated fragments would yield some 
correlation. In other words, projected HFB is not size 
consistent since it does not dissociate to projected HFB 
fragments. In the special case of NRHFB = AGP, this 
is well known.— The lack of size consistency is not fa- 
tal for molecular applications, but is of course a serious 
problem for applications to periodic systems, where the 
correlation energy per unit cell would vanish as the size 
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FIG. 7. Dissociation of He2 in the 3-21G basis. We have 
broken the energy axis several times as the solutions are well- 
separated in energy and the binding is minimal. The full CI 
results exactly match the KNRHFB curve shown in the figure. 
The labels "s" and "bs" on the NRHFB indicate whether 
spatial symmetry is preserved ("s") or broken ("bs") in the 
underlying HFB determinant. 



of the unit cell is increased. 

We consider the dissociation of He2 in Fig. [71 using the 
3-21G basis. The RHF curve is strictly repulsive, while 
NRHFB has noticeable binding strictly because it goes 
to the wrong dissociation limit. This binding is greatly 
exaggerated, as the experimental binding energy is on 
the order of 30 microHartrees. By breaking spatial sym- 
metry, we obtain a curve that closely parallels the RHF 
curve; in fact, the dissociation limit corresponds to one 
RHF atom and one AGP atom. This suggests that spa- 
tial symmetry restoration (here equivalent to restoring 
inversion symmetry due to the small basis being used) 
has something to offer. In fact, when we restore spa- 
tial symmetry in Ci-NRHFB, we obtain a curve which 
is only about 0.1 milliHartree above full CI. Curiously, 
restoring instead complex conjugation symmetry in KN- 
RHFB gives results which are numerically identical to 
the FCI. 

Comparing NRHFB to C; NRHFB, we note that one re- 
covers all but ~ 0.1 milliHartree of the ~ 15 milliHartrees 
error in the dissociation limit of the broken symmetry 
NRHFB. reducing the error by more than 99%. We sus- 
pect that by breaking and restoring the remaining possi- 
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FIG. 8. Correlation energies in the four-electron series. We 
compare the NRHFB results in the aug-cc-pVTZ basis to the 
exact results of Ref. [stI . Note that we correctly recover the 
linear behavior of the correlation energy as a function of nu- 
clear charge Z. 



ble symmetries we may be able to achieve FCI results. 



F. Accounting for Static Correlation 

In previous plots, we showed that the projected HFB 
wave function can describe left-right correlation, a form 
of strong correlation that occurs in molecular dissocia- 
tion to fragments. This can be accomplished by break- 
ing the spatial and spin symmetry on top of the number 
projected wave function (as in the case of the NUHFB 
dissociation of N2), or by breaking and restoring other 
symmetries (as in the case of the KNRHFB dissociation 
of H4). 

Another type of strong correlation is referred to as 
angularj^ and is related to the near degeneracies of or- 
bitals with the same principal quantum number in atomic 
systems. For instance, the 2s and the 2p orbitals of beryl- 
lium are close in energy, and thus a multireference wave 
function is needed to give a correct description of the 
electronic structure of this atom. 

We show in Fig. [8] the correlation energy predicted 
by the NRHFB method on the series of atomic species 
having 4 electrons. We compare to the exact correla- 
tion energies at the non-relativistic, Born-Oppenheimer 
level of theory provided by Chakravorty et al.—. We 
have used the aug-cc-pVTZ basis for NRHFB calcula- 
tions; this basis is large enough to provide qualitatively 
correct behavior, but still departs from the complete ba- 
sis set limit. The leading contribution to the correlation 
energy is known to scale linearly with Z for large enough 
and NRHFB reproduces this behavior. 

We note that if the core occupations were frozen (both 
in FCI and in NRHFB), then NRHFB would agree ex- 



actly with FCI because what remains is a two-electron 
system and AGP (or NRHFB) is exact for any two- 
electron singlet. 

G. Accuracy of Projections 

In order to assess the number of grid points required in 
the particle number projection, we have performed cal- 
culations on rings of hydrogen atoms of increasing size. 
In our results, we have converged NRHFB states with 9 
grid points in the integration from to tt using the trape- 
zoidal rule and the equations derived in Ref. 0. (Since 
the number of particles is even, the integrand is sym- 
metric about 77.) We then evaluated errors in the energy 
and in the number of particles with grids of various sizes 
using the converged density matrix. 

Our results are shown in Table [Hi We note that the 
trapezoidal rule becomes exact when the number of grid 
points used satisfies M — max{^N, fl— ^N) + l, where 
is the total number of particles and D, is the degeneracy 
of the spacci^ Thus, the trapezoidal grid becomes exact 
for Hg in minimal basis with 5 grid points. 

The numbers shown in Table |ll] suggest that the size of 
the grid required to achieve certain accuracy in the en- 
ergy or in the variance of the number of particles scales 
less than linearly with the number of particles. We expect 
the same weak dependence for other projection operators. 
Presumably, the quality of the grid required for spin pro- 
jection depends on the magnitude of the spin contamina- 
tion present in the broken symmetry HFB state. 

V. CONCLUDING REMARKS 

Symmetry breaking plays a central role in our cur- 
rent understanding of classical and quantum systems.— 
In finite electronic structure systems, spontaneous sym- 
metry breaking in the presence of exact or near degen- 
eracies is an artifact and quantum fluctuations (strong 
correlations) remove them. In the presence of degen- 
eracies, strong correlation therefore cannot be neglected 
if one wants to obtain a correct qualitative picturei^ 
One should remark that Schrodinger's equation is linear 
whereas the mean field equations of HF and HFB theories 
are cubic in the orbitals, and have more solutions than 
the physical ones. Spontaneous symmetry breaking in 
mean-field theories flags emerging behavior, the appear- 
ance of phenomena (due to degeneracies) that requires 
a description beyond the single determinant picture. By 
breaking symmetry, mean-field theories lift the degen- 
eracy {e.g. the HOMO-LUMO gap opens going from 
RHF to UHF in H2 near dissociation) and can some- 
times give qualitative descriptions of these phenomena at 
the cost of good quantum numbers. In this sense, mean- 
field theories predict their own failure and signal the need 
for a more comprehensive treatment. Projection after 
variation is not the answer because orbital relaxation 
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TABLE II. Errors associated with the discretization of the particle number projection operator in NRHFB calculations of 
hydrogen atom rings (rn-H ~ 1-80 bohr, ST0-3G basis). The reference energy is taken to be that with A'^grid = 9. 



Property 


AT 










abs ) — iv j 


I 


/.OO X iU 


6.U0 X iU 


l.ZD X iU 


0.(1 X iU 




4 


4.88 X 10 


1.46 X 10 


1.43 X 10 


9.06 X 10 




6 


< 10-1° 


1.41 X IQ-'' 


6.95 X IQ-^ 


9.73 X 10-^ 




8 


< 10-1° 


< 10-1° 


< 10-1° 


< 10-1° 


abs [(^Jy ) — [Jy) ) 


2 


2.ZI X 10 


Z.oZ X iu 


o ^ / in — 1 
2.87 X iu 


o OT V/ in — 1 

z.o7 X 10 




4 


7.03 X 10-^ 


3.55 X 10"* 


7.99 X IQ— * 


1.17 X 10-^ 




6 


< 10-1° 


2.41 X IQ-* 


2.20 X IQ-'' 


6.90 X lO-'^ 




8 


< 10-1° 


< 10-1° 


< 10-1° 


1.15 X 10-1° 


ahs{{H) - E,,t) 


2 


5.92 X 10-^ 


5.46 X IQ-^ 


4.78 X 10-^ 


4.13 X 10-^ 




4 


1.54 X 10-^ 


5.91 X IQ-^ 


9.69 X 10-^ 


1.10 X 10-1 




6 


< 10-1° 


3.53 X IQ-'' 


2.32 X 10-* 


5.33 X 10-** 




8 


< 10-1° 


< 10-1° 


< 10-1° 


< 10-1° 



effects are very important. With PAV, the unphysical 
mean-field behavior is frequently enhanced rather than 
ehminatedi ^^i'^^ On the other hand, the results in this 
paper seem to indicate that a comprehensive variation- 
after-projection treatment of all molecular symmetries is 
capable of accounting for molecular static correlation in a 
black-box manner that yields smooth dissociation curves. 
The fact that this can be achieved with mean field compu- 
tational cost is truly remarkable and unprecedented. On 
the other hand, extended systems behave differently from 
finite systems and it is usually argued that spontaneous 
symmetry breaking there has the physically meaningful 
interpretation associated with true phase transitions.— 

The theory presented in this paper covers important 
unexplored aspects of electronic structure theory. The 
first one is the underlying coherent state representation 
where the generator coordinate method projects out vari- 
ational states of the correct symmetry>ii^ The manifold 
of states from which this projection is performed is non- 
orthogonal and overcomplete. Our projected quasipar- 
ticle states are multirefcrence wave functions. A second 
salient aspect of our work is the connection with gcm- 
inal theories 1^ We have extended these geminal wave 
functions to variationally include unrestricted and gen- 
eral spin orbitals. We have also calculated wave func- 
tions that are linear combinations of these general-orbital 
AGPs. 

Compared to CPMFT, our previous model for strong 
correlations, the current theory is N-representable and 
possesses a two-particle density matrix that is factoriz- 
able over the gauge grid (see Appendix C). As a matter of 
fact, all higher order reduced density matrices factorize in 
a similar way. This property will certainly be of interest 
to other workers who are building correlation models on 
top of multirefcrence descriptions that do not possess this 
very useful factorization, since the factorization leads to 



one power reduction in computational scaling (typically 
from 0(M6) to 0{M^)). 

Regarding the separation of static and dynamic corre- 
lation that we have previously advocated in our CPMFT 
workjii-— one should note that the present theory is ex- 
act for any two-electron system. This, in our previous 
definition, includes both static and dynamic correlation. 
In the present geminal context, it seems more advanta- 
geous to describe correlations as originating from inter- 
and intra-pair interactions. It is evident that PHFB in- 
cludes some but not all dynamical correlations although 
it seems to include all strong correlations for molecules. 

One more time, recapitulating the main points of this 
work: 

• First, we work with underlying unprojected quasi- 
particle determinants which deliberately break 
symmetries of the exact wave function. These 
symmetries are restored by projection to give pro- 
jected quasiparticle wave functions which are ob- 
tained variationally. We are following the proce- 
dure, that is, of variation-after-projection. The 
projected quasiparticle wave functions are multiref- 
crence in character. 

• Second, because the underlying unprojected wave 
function is simply a single determinant of quasi- 
particles, and the projected wave function is com- 
pletely specified by the unprojected wave function, 
our projected wave function is specified by a regular 
density matrix p and an anomalous density matrix 
K. The problem of variation after projection re- 
duces to the problem of optimizing p and k for the 
unprojected state, which can be accomplished at 
mean-field computational cost. 

• The projection operators we are using are written 
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as integrations over gauge angles in a generator co- 
ordinate approach. We simply discretize these in- 
tegrals to obtain numerically efficient projections, 
and fortunately the grids we need for each projec- 
tion operator are not large. We should point out, 
also, that while the computational cost is roughly 
equivalent to that of a mean field calculation at ev- 
ery grid point, the problem of gauge integration to 
do the projections is trivially parallel. 

In summary, to the best of our knowledge the main 
new accomplishments presented here are: 

• We have presented an efficient algorithm allowing 
for AGP calculations at mean field computational 
cost. 

• We have presented the first variational AGP calcu- 
lations with broken symmetry orbitals. 

• We have reported the first full VAP calculations 
based on HFB for the combination of number, spin, 
and complex conjugation restoration. 

• We have reported the first discrete symmetry 
restoration (complex conjugation and point group) 
in quantum chemistry. 

• We have reported the first application of VAP with 
the full electronic Hamiltonian. 

In closing, we would like to emphasize that the calcu- 
lations presented in this work are just proof-of-principle 
benchmarks and only meant to demonstrate the com- 
pelling capability of the theory. Larger bases and chem- 
ically meaningful results will be presented in due time. 
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Appendix A: Equations for the Effective 
Hamiltonian 



Here we present the expressions for the effective Hamil- 
tonian in the case of a general projection operator P as- 
sociated with the rotation operator R{9). We note that 
in their 2000 paper. Sheikh and Ring derived a simpli- 
fied form of these equations applicable only to particle- 
number projection. 

Let us begin by discussing the general form of the ef- 
fective Hamiltonian. Recall that the energy is given by 



E^IJ d9y{9)Tr[{h + Te)pe 



which we will express equivalently as 

E^IJ dey{0)TT[{2h + Ge)pe- ^e^*s]. (A2) 
The effective Hamiltonian will take the general form 



pp + A'' + A*" 

-(F^ + A-")* -{FP + AP)* 



where 



pp. 



d 1 



dpji 
d 1 



A". 



A 



dp,, 2 
d 1 



\ J ddy{d)Tii{2h + Ge)pel 
-IJ d0y{e)Tr[{2h + Gg)pg] 
d0y{9)Tv[AeK*], 



d 1 



d0yid)TT[AeK*g] 



(A3) 

(A4a) 
(A4b) 
(A4c) 
(A4d) 



Dependence on p and k lurks in y{9) and in the transition 
density matrices pe, ttg, and Kg (and thus in Gg , Ag , and 

Ag). 

Following Sheikh and Ring, we write the derivatives of 
the function y{6) as 



dyjO) 
dpki 
dyjO) 



= y{e)Yii{e) = y{e)[Yik{e)-%i{e)]. 



pp. 



We therefore obtain 

dey{e)[YPTv[{h + \Gg)pg\ 

dey{e)[Y:^TT[{h+]^Gg)pg\ 



A". 



\ j dey{e)[YPTv[AgKl] 



A..(^)^ + Ar,(.)^ 



dpjt 

^ J d0yie)[Y{;TiiAg^*, 



(A5) 
(A6) 

(A7a) 
(A7b) 
(A7c) 

} 

(A7d) 

} 



(Al) 



We emphasize that PP, F'^, Ap, and A'' are general- 
izations of the usual Fock and pairing matrices of HFB 
and CPMFT. The derivatives of the transition density 
matrices are involved and we merely provide the final 
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expressions. Note that F^, A'', and V are Hermitian, 
though this is not obvious from the matrix expressions 
below and we therefore Hermitize at the end. Similarly, 
pK ^ A'', and V are antisymmetric and we have explic- 
itly antisymmetrized them in the course of taking the 
derivatives. 

Our final results are that 

= i (RepR^Ce + RjCepRe) (A8) 

- i y <:\^y{4>) (R^pR^C^ + R^C^pR^) 

Ye = -iRtCe/tR? ^\ j d0y((/))R;C^/^R; (A9) 
Y« = Ye - Ye^ (AlO) 
Using these, we have 



pp 



\\ Ue y{e) \yp TT[{h + ]^Gg)pg] (Alia) 



(1 - pe):^eRepR^Ce 



h.c. 



dey{6) YTr[{h+-Ge)pe] (Allb) 



Evaluating the numerator requires integration over the 
projection operator grid twice. We have used this ap- 
proach, for example, in evaluating density matrices (see 
below). For the special case that the projection opera- 
tor commutes with the operator O, the expectation value 
simplifies to 



(0) = 



($|0P|$) 



(B2) 



as we have used to evaluate the energy. 

Because the spin projection operator takes the form 
exp(iaiS'z) exp(i/3S'y) exp(i7iS'z), it does not commute 
with the individual Si, so evaluating their expectation 
values is most conveniently done in terms of the one- 
particle density matrix given below. However, since Sy 

and 5*2 both commute with , evaluation of {S"^) is sim- 
pler. We find that 

(5^)= j d02;(0)|^(^Tr[M,(0)]2 (B3) 



+ -Tr[M2(0)-/^*(0)/^, 
+ ^Tr[P(0)-p2(0)-KS(0)Ko( 



where we have decomposed the spin blocks of the transi- 
tion density matrices as 



A" 



d^^2/(0)[- Y^TriAe/^^] 

R^CepAeKgRe 

(1 - pe)^gIiWB)gC0 

RtCeKAJ(l - pe)Re 

KffAeRepRjCfl ^ + h.c. 



dey{e) [YTr[Ae/^^] 
RjCepAepjR; 
RtCeKAgKeRfl 



(Allc) 



(Alld) 



Appendix B: Evaluation of (S^) 

In general, the expectation value of an operator O with 
our PQT wave functions is 



(0) = 



($|POF|$) 



(Bl) 



P{0) 



(B4) 



<0) 



(B5) 



,ppa{o) ppm) 

' P((?)+M,((?) M,(0)-iM,(0)^ 
^M,(0)+iM,(0) P(0)-M,(0) ^ 



Appendix C: Form of the Reduced Density Matrices 

To compute a general expectation value, it is most 
convenient to work with reduced density matrices. We 
discuss here how to evaluate these tensors using the pro- 
jected HFB state. The A^-particle reduced density matrix 
(A^-RDM) corresponding to the projected HFB state is 
given by 



Nr ^1 



■■lN,Jl---]N 



($|ptp|$) 



(Cl) 
(C2) 
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where we have defined 



with 



^il---iN ,jl---jN ~ jy-j ^ji ' ' ' '^jjv'^'N ' ' ' '^^i ■ (^3) 

Unfortunately, the operators F do not commute with 
the projection operators in general. However, one can 
evaluate the reduced density matrices by double integra- 
tion over the gauge angle. That is, 

P ^ nded9'w{e)w{e')mRH9')tR{9m 
jjd0de'w{e)w{e')mRf{e')R{9)\<P) 

Defining the normalized weighting function 

// d(/) d(/.' ) I (0) I 

the reduced density matrix becomes 



dode' yie,e'){e'\r\e) 



(C6) 



The overlap matrix elements can be computed from 

det Rfl det Rt 



{^\R\9') R{9)\^) = 



^/det p-y/det Cg 



(C7) 



Cgg, — Rg' p Rg, Re p Rg — Re' k RJ, Rg k* Rg. (C8) 

The transition density matrices needed to evaluate the 
density operator expectation values are in turn given by 

pi9,9') ^RepRlCee'Ro' pRln (C9) 

K{9,9')=RgpRlCgg,Rg,KRl, (CIO) 

K*{9,9') ^R*gK*RlCge'Rg' pRl,. (CU) 

The iV-PDM is then formed by integrating products of 
transition density matrices. For example, the 2-PDM is 



Tkl,^J=l II d9d9' 



'yi9,9')[pk^{9,9')pi,i9,9') (C12) 
-Pkj{9,9')pu{9,9') 
+ nt^{9,9')Kki{9,9')} 



with similar factorizable expressions for all higher order 
density matrices. 
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